In this lab you will learn



library(tidyverse)
library(gridExtra) ## for plotting function grid.arrange()
library(magrittr) ## for pipe symbol %<>%
library(knitr) ## for table-viewing function kable()
library(matrixcalc) ## for function matrix.power()
library(rio) ## for function import to load Excel dataset from GitHub
theme_set(theme_bw())
theme_update(panel.grid = element_blank())



All living organisms go through life stages as they develop. In most species, those life stages differ in ecologically important ways. Most notably, both the probability an individual survives until a later stage and their ability to produce offspring can change considerably over the course of their lives.

Consider for example the table of vital statistics below for American beech, Fagus grandifloia.

Seedlings of beech and maple in a beech–maple forest, by Quercus2018


A simplified tree life table for American beech in southeast Texas (data from Harcombe and Marks 1983 via Harcombe 1987).
Size class Annual proportion dying Annual proportion growing into the next size class Annual per capita seed production
Seeds 0.90 0.10 0
Seedlings (<50 cm tall) 0.65 0.05 0
Saplings (<50 cm tall to 4 cm dbh 0.08 0.02 0
Poles (4-30 cm dbh) 0.06 0.01 0
Mature trees (>30 cm dbh) 0.02 NA 200


Notice how the probability a seedling survives to sapling stage is only 35%, while 94% of the poles survive to become mature trees. In contrast, only mature trees reproduce, while all other stages have zero fertility!

As another example, consider human fertility by age of the mother in the United States:


Table: Total fertility rates are the sums of birth rates by 5-year age groups of mothers multiplied by 5. Birth rates are live births per 1,000 women in specified group. Data for 2003 from the CDC.

This marked difference in vital rates – fecundity and mortality – across life stages has important consequences for the growth and eventual makeup of the population. In this lab we will explore how differences in vital rates across ages or stages as in the examples above affect the age- or stage-structure in the population.


Age-structured growth

As we have seen in lecture, age-structured populations can be described by a projection matrix with a very specific structure, because surviving individuals always move to the next age class between time steps, and new individuals can only enter the population in the first age class. This matrix is named the Leslie matrix after the British ecologist Patrick H. Leslie. It contains mostly 0’s, except for the fecundity rates (number of offspring per capita per time step) and the survival rates (proportion of individuals who survive and move to the next age class) of each age class.

For a population with \(k\) age classes, the Leslie matrix takes the form

\[ \begin{pmatrix} f_1 & f_2 & f_3 & \dots & f_k\\ s_1 & 0 & 0 & \dots & 0\\ 0 & s_2 & 0 & \dots & 0\\ 0 & 0 & s_3 & \dots & 0\\ \vdots & \vdots & \vdots & \ddots & 0 \\ 0 & 0 & 0 & s_{k-1} & 0 \end{pmatrix} \] where \(f_i\) and \(s_i\) are respectively the fecundity and survival rate of age class \(i\). In words, the only possibly non-zero entries are the first row and those immediately below the diagonal.


Cycling in age-structured populations

Consider a population of beetles with two age classes: juveniles (larvae + pupae) and adults. The species spends the first year of its life as a juvenile, and then those who did not die from predation or disease become adults. The average adult beetle of that species produces 10 hatching eggs, and all adults die within a year of emerging from their cocoons. You estimate from field observations that this species suffers 90% mortality at the juvenile stage.

The Leslie matrix for this population, since it has two classes, will have two rows and two columns: the first row (column) corresponds to the juveniles, and the second to adults. The matrix will look like this:

\[ \begin{pmatrix} a & b \\ c & d \end{pmatrix} \]

Q1: Based on the information about the beetles given above, what should be the values of \(a\), \(b\), \(c\), and \(d\)?


Q2. Consider an initial beetle population of 100 juveniles and 20 adults. Modify the code below in the appropriate places to specify the Leslie matrix and initial conditions, then run a simulation of the dynamics of this structured population and show a plot of the ratio of juveniles to adults over time. Describe what you see.

## Leslie matrix
L = matrix(c(a, b, c, d), byrow = TRUE, nrow = 2)

## initial conditions
N = c(J0, A0)
N_record = N

## number of years to run for
years = 10

## dynamics
for(year in 1:years){
  N = L %*% N
  N_record = cbind(N_record, N)
} 

# collect results in a tibble
N_tibble = 
  tibble(
    year = 0:years, 
    juveniles = N_record[1, ],
    adults = N_record[2, ],
    ratio = juveniles / adults
  )

#plot
N_tibble |> 
  ggplot(aes(year, ratio)) + 
  geom_point() +
  geom_line()


Let us try to understand the behavior of the beetle population. The Leslie matrix in the population above is of the form

\[ \pmatrix{ 0 & f \\ s & 0 } \]

where \(m\) is the fecundity of adults and \(p\) is the survival of juveniles (recall that adults have zero survival and juveniles have zero fecundity). Suppose in Year 0 the population has \(J_0\) juveniles and \(A_0\) adults, which we can write in matrix notation as the state

\[ \bf{N}_0 = \pmatrix{J_0 \\ A_0} \]

The next year, the population will be

\[ \bf{N}_1 = \pmatrix{ 0 & f \\ s & 0 }\pmatrix{J_0 \\ A_0} \]

In Year 2, the population will then be

\[ \bf{N}_2 = \pmatrix{ 0 & f \\ s & 0 }\pmatrix{ 0 & f \\ s & 0 }\pmatrix{J_0 \\ A_0} \]

We can simplify this expression by working out the matrix product:

\[ \pmatrix{ 0 & f \\ s & 0 }\pmatrix{ 0 & f \\ s & 0 } = \pmatrix{sf & 0 \\ 0 & sf} = sf \pmatrix{1 & 0 \\ 0 & 1} \]

This means that

\[ \bf{N_2} = \it{sf} \pmatrix{1 & 0 \\ 0 & 1} \pmatrix{J_0 \\ A_0} \]

You may recognize the matrix \(\pmatrix{1&0\\0&1}\) as the identity matrix, which is essentially the analogue of the number 1 in matrix world: just as 1 times any number \(n\) equals \(n\), the product of the identity matrix with any vector gives the same vector. We conclude that

\[ \bf{N}_2 = \it{sf} \pmatrix{J_0 \\ A_0} \] In other words, in Year 2 there will be \(J_2 = sf \, J_0\) juveniles and \(A_2 = sf \, A_0\) adults.


Q3. What will be the ratio of juveniles to adults in Year 2? How does that compare to the ratio of juveniles to adults in Year 0?

Projecting this forward, we now can see that in Year 3 the population will look similar to Year 1, then in Year 4 it will again look like Year 0, and so on. Which is to say, the demographic structure of this population is locked in a two-year cycle.

What causes the cycling? We can figure it out by thinking about what’s happening: Every year, a proportion \(p\) of the previous cohort of juveniles becomes the new cohort of adults. One year later, each of those adults produces \(m\) new juveniles. But this means that the size of this new cohort of juveniles is proportional to the size of the cohort of juveniles from two years ago, and the constant of proportionality is \(sf\). You can apply the exact same rationale to convince yourself that the same is true of what’s happening to the adult cohorts.

Note: Mathematically, this all happened because the Leslie matrix of that population had a special property: when squared, it becomes similar to the identity matrix. In general, the Leslie matrix of a population with \(k\) stages where only the last stage reproduces will become similar to the identity matrix after being multiplied by itself \(k\) times, leading to a cycle every \(k\) time units. We saw an example of this periodic cycling in lecture.

Anytime we make an intriguing theoretical discovery, we must ask ourselves whether the result is robust to dropping unrealistic assumptions. Otherwise, our theoretical finding will have limited real-world implications. In our case, we should check that the cycling will happen even if we allow adult fecundity and juvenile survival to vary between the years, as it is not realistic to assume they will be fixed.


Q4. The code below repeats the experiment from Q2, but this time we allow adult fecundity and juvenile survival to fluctuate from one year to the next (this is called environmental stochasticity in ecological jargon). Run the code, then plot the ratio of juveniles to adults by running the command plot_ratio and describe the results. Does the ratio still cycle? What is different from the fixed-rates scenario above? (If you are wondering if the observed behavior eventually stops, increase the number of years that we run the dynamics for and check out the plots.)

# control random number generator
set.seed(0)

# vital rates
time_avg_m = 10
time_avg_p = .1


# initial conditions
J0 = 100
A0 = 20
N = c(J0, A0)
N_record = N

# number of years to run for
years = 10

# dynamics
for(year in 1:years){
  
  ## Leslie matrix
  m = max(0, rnorm(1, mean = time_avg_m, sd = .1 * time_avg_m))
  p = max(0, rnorm(1, mean = time_avg_p, sd = .1 * time_avg_p))
  L = matrix(c(0, m, p, 0), byrow = TRUE, nrow = 2)

  N = L %*% N
  N_record = cbind(N_record, N)
} 

# collect results in a tibble
N_tibble = 
  tibble(
    year = 0:years, 
    juveniles = N_record[1, ],
    adults = N_record[2, ],
    ratio = juveniles / adults
  )

# plot
plot_ratio = 
  N_tibble |> 
  ggplot(aes(year, ratio)) + 
  geom_line() +
  geom_point() +
  ylab('Ratio of juvenile population to adult population')


It is not just the ratio of juveniles to adults that cycles. The plot below shows the adult population from the simulation above.

N_tibble |>
  filter(year < 10) |>
  ggplot(aes(year, adults)) +
  geom_line() +
  geom_point(size = 3) +
  theme(aspect.ratio = 1)


Q5. Show the analogous plot for juveniles.

The results above suggest that we should expect some age-structured populations to cycle in the real world. Indeed, population cycling is a well-documented phenomenon. The figure below was taken from a study on the red flour beetle, Tribolium castaneum. It shows clear cycling in the populations of larvae, pupae, and adults.

Triboleum castaneum, the red flour beetle

Time-series data (black dots) on experimental populations of T. castaneum show cycling at the larval (L), pupal (Projection_matrix), and adult (A) stages. The unit of time is two weeks. The constant \(c_{pa}\) on the plots refers to the adult-pupa cannibalism coefficient. From Figure 9 in Dennis et al. 2001, Ecological Monographs 71 (2) 277-303.



Population pyramids

As we have seen in lecture, we can use the Leslie matrix to predict the stable age distribution of a population. Here we will make such prediction for women in the United States.

Our Leslie matrix will contain the fecundity rates and survival probabilities of 5-year age groups. We will use survival data for 2008 sourced from the CDC/NCHS National Vital Statistics System, and fecundity data for 2003 sourced from the CDC.

github_address = 'https://github.com/rafaeldandrea/BIO-356/blob/master/Survival_table_CDC.xlsx?raw=true'

survival = 
  rio::import(github_address, range = 'A3:D25') |> 
  transmute(age = `...1`, survivors = Female) |>
  filter(age != 1) |>
  as_tibble() |>
  mutate(survival = dplyr::lead(survivors) / survivors)

fertility = 
  tibble(
    `Age of mother` = 
      c('10-14', '15-19', '20-24', '25-29', '30-34', '35-39', '40-44', '45-49'),
    Fertility = c(.6, 41.6, 102.6, 115.6, 95.1, 43.8, 8.7, .5)

  )

vital_usa =
  survival |>
  left_join(
    fertility |>
      transmute(
        age = seq(10, 45, by = 5),
        fertility_woman = Fertility,
        fertility_group = fertility_woman / 1000 * 5 * .4878 
      ),
    by = 'age'
  ) |>
  replace_na(list(fertility_woman = 0, fertility_group = 0))


Table: Age-specific vital statistics for women in the United States. Age group: Individuals within five years of the listed age. Eg, age group 30 contains all women aged 30-34. Survivors: Number of survivors out of 100,000 born alive. Survival: proportion of individuals in an age group who survived into the next age group, ie survival(age group) = survivors(next age group) / survivors(age group). Fertility per 1,000 women: number of live births of both sexes per 1,000 women in the age group. Fertility (daughters): number of female live births per woman in age group, obtained from the previous column F via the formula m = F / 1000 * 0.4878, the latter constant being the female-to-male sex ratio at birth in humans.
Age group Survivors Survival Fertility per 1,000 women Fertility (daughters)
0 100000 0.993 0.0 0.000
5 99305 0.999 0.0 0.000
10 99249 0.999 0.6 0.001
15 99185 0.998 41.6 0.101
20 99022 0.998 102.6 0.250
25 98794 0.997 115.6 0.282
30 98519 0.996 95.1 0.232
35 98169 0.995 43.8 0.107
40 97665 0.992 8.7 0.021
45 96857 0.987 0.5 0.001
50 95606 0.981 0.0 0.000
55 93810 0.973 0.0 0.000
60 91317 0.959 0.0 0.000
65 87571 0.937 0.0 0.000
70 82039 0.902 0.0 0.000
75 73974 0.844 0.0 0.000
80 62448 0.749 0.0 0.000
85 46782 0.594 0.0 0.000
90 27805 0.396 0.0 0.000
95 11006 0.213 0.0 0.000
100 2345 NA 0.0 0.000


The \(s\) entries in our Leslie matrix come from the third column in the table above (Survival), and the \(f\) entries come from the last column above (Fertility (daughters)).

As discussed in lecture, once we have the Leslie matrix, we can use its leading eigenvector to predict the stable age distribution.

The code chunk below defines the function Predicted_Population_Pyramid, which takes as input the vital statistics and total population size, and outputs a plot of the predicted population pyramid.

Note: The following is a bit technical but helps you understand the code: In R, we can very easily retrieve the eigenvectors of a matrix using the function eigen(), which also returns the eigenvalues. We extract the eigenvectors using the identifier $vectors and then single out the leading eigenvector with the specifier [, 1]. Finally, we discard the imaginary part of the leading eigenvector (which is zero anyway) using the function Re(), and then wrap it in the function abs() to make sure our eigenvector components are positive. The last step is needed because R returns them as either all positive or all negative, since multiplying all components by -1 (or any other constant) does not affect the eigvenvector.


Predicted_Population_Pyramid = function(vital_statistics, total_population){
  # create matrix filled with 0's
  Leslie_matrix = matrix(0, nrow = nrow(vital_statistics), ncol = nrow(vital_statistics))
  
  # fill in the fecundity data on the first row
  Leslie_matrix[1, ] = vital_statistics$fertility_group
  
  # fill in the survival data below the diagonal
  for(j in 1:(nrow(vital_statistics) - 1)) Leslie_matrix[j + 1, j] = vital_statistics$survival[j]
  
  # calculate eigenvalues and eigenvectors
  eigenvalues_and_eigenvectors = eigen(Leslie_matrix)
  
  # extract eigenvectors 
  eigenvectors = eigenvalues_and_eigenvectors$vectors
  
  # subset leading eigenvector
  leading_eigenvector = eigenvectors[, 1]
  
  # take positive real part
  leading_eigenvector_real_and_positive = abs(Re(leading_eigenvector))
  
  # renormalize eigenvector to sum = 1
  scaled_eigenvector = 
    leading_eigenvector_real_and_positive / sum(leading_eigenvector_real_and_positive)
  
  # estimate the predicted stable age distribution as the components of the leading eigenvector,
  # normalized by the total female population in the US in 2020 (.508 of 328.2m people).
  predicted_pyramid =
    tibble(
      age = factor(vital_statistics$age), 
      predicted_population = scaled_eigenvector * total_population
    )
  
  # plot the predicted population pyramid
  plot_pyramid =
    predicted_pyramid |>
    ggplot(aes(age, predicted_population)) +
    geom_col(fill = c('salmon', "#DD5144")[1 + seq(nrow(Leslie_matrix)) %% 2]) +
    coord_flip() +
    theme(aspect.ratio = 1) +
    labs(
      y = 'Population (millions)',
      x = 'Age group'
    )
  
  plot_pyramid
}


Now we are ready to plot our predicted stable age distribution by calling the function Predicted_Population_Pyramid defined above, with arguments vital_statistics = vital, and total population = 331 * 0.508, which is the estimated 2020 female population of the Unites States.

Predicted_Population_Pyramid(vital_statistics = vital_usa, total_population = 331 * .508)

We can compare our predictions with the actual population pyramid of the United States:

Population pyramid of the US, by CIA Factbook

Not a bad match considering the simplicity of our age-structured model!

Let’s now consider another country, Nigeria. The fertility and survival curves in Nigeria differ substantially from their US counterparts, as shown below.

Data: Fertility from Nigeria Demographic and Health Survey 2018 Table 5.1, survival from World Health Organization.

survival_nigeria =
  tibble(
    age = seq(0, 85, by = 5),
    survivors = 
      c(
        100000,
        89237.13,
        88086.45,
        87489.52,
        86035.28,
        83883.13,
        81529.16,
        78674.88,
        75099.65,
        70980.3,
        66487.29,
        61201.56,
        54937.09,
        46545.7,
        36142.82,
        23885.7,
        12046.73,
        3965.038
      )
  ) |>
  mutate(survival = dplyr::lead(survivors) / survivors)

fertility_nigeria = 
  tibble(
    age = seq(10, 45, by = 5),
    fertility_woman = c(2, 106, 239, 256, 217, 149, 67, 23),
    fertility_group = fertility_woman / 1000 * 5 * .4878
  )

vital_nigeria = 
  left_join(
    survival_nigeria,
    fertility_nigeria,
    by = 'age'
  ) |>
  replace_na(list(fertility_woman = 0, fertility_group = 0))


Table: Age-specific vital statistics for women in Nigeria. Columns are analogous to those in the USA table above.
Age group Survivors Survival Fertility per 1,000 women Fertility (daughters)
0 100000.000 0.892 0 0.000
5 89237.130 0.987 0 0.000
10 88086.450 0.993 2 0.005
15 87489.520 0.983 106 0.259
20 86035.280 0.975 239 0.583
25 83883.130 0.972 256 0.624
30 81529.160 0.965 217 0.529
35 78674.880 0.955 149 0.363
40 75099.650 0.945 67 0.163
45 70980.300 0.937 23 0.056
50 66487.290 0.921 0 0.000
55 61201.560 0.898 0 0.000
60 54937.090 0.847 0 0.000
65 46545.700 0.777 0 0.000
70 36142.820 0.661 0 0.000
75 23885.700 0.504 0 0.000
80 12046.730 0.329 0 0.000
85 3965.038 NA 0 0.000


Q6. a) Use the code below to plot the survival and fecundity curves for women in the US and Nigeria. b) For each vital rate, describe the differences and similarities between the two countries. c) Given the differences, how should the Nigerian population pyramid differ from that of the US? (Note: To run the code below you must first run the code chunks reading the data for the US and Nigeria.)

vital_data = 
  bind_rows(
    vital_usa |> mutate(country = 'USA'),
    vital_nigeria |> mutate(country = 'Nigeria')
  )

plot_survival = 
  vital_data |>
  filter(age <= 80) |>
  ggplot(aes(age, survival, group = country, color = country)) +
  geom_line() +
  geom_point() +
  theme(
    aspect.ratio = 1,
    legend.position = 'top'
  ) +
  labs(
    x = 'Age group',
    y = 'Survival'
  )

plot_fecundity = 
  vital_data |>
  filter(age <= 80) |>
  ggplot(aes(age, fertility_group, group = country, color = country)) +
  geom_line() +
  geom_point() +
  theme(
    aspect.ratio = 1,
    legend.position = 'top'
  ) +
  labs(
    x = 'Age group',
    y = 'Live female births per woman in age group'
  )

gridExtra::grid.arrange(plot_survival, plot_fecundity, nrow = 1)


Q7. Call the function Predicted_Population_Pyramid using the Nigerian demographic data to plot the population pyramid of women in Nigeria (total number of women is approximately 99.13 million). Describe the differences from the US counterpart. If the result differs from your prediction in Q6, connect those differences to the differences in the fecundity and survival curves of the two countries. (Note: Before you can call the function, you must first run the code defining the function. Find and copy/paste the relevant code chunk above)



Stage-structured populations

Consider the life table of American beech shown at the top of this lab write-up. This is a stage- or class-structured population, not an age-structured population, as members of each class may remain in the same class between time units. They could also in principle regress to earlier stages, although that does not happen in this particular population.

Figure: Fagus grandifolia fall leaves, by Katja Schulz, NCSU

Based on the fecundity and mortality rates of the different life stages, we can write its projection matrix as follows:

Projection_matrix = 
  matrix(
    c(
       0,   0,   0,   0, 200, 
      .1,  .3,   0,   0,   0, 
       0, .05,  .9,   0,   0,
       0,   0, .02, .93,   0,
       0,   0,   0, .01, .98  
    ),
    byrow = TRUE,
    nrow = 5
  )

Projection_matrix
##      [,1] [,2] [,3] [,4]   [,5]
## [1,]  0.0 0.00 0.00 0.00 200.00
## [2,]  0.1 0.30 0.00 0.00   0.00
## [3,]  0.0 0.05 0.90 0.00   0.00
## [4,]  0.0 0.00 0.02 0.93   0.00
## [5,]  0.0 0.00 0.00 0.01   0.98


Q8. Assuming unchanging vital rates, is this population projected to grow or decline in the long term? To find out, run the code chunk above for the projection matrix, then use the command line leading_eigenvalue = Re(eigen(Projection_matrix)$values[1]).


Q9. Based on the entries in the projection matrix, which stages do you predict to be most and least abundant in the long run? Check how your intuition compares against the predicted stable class distribution using the code below. Plot the stable class distribution.

right_eigenvector = abs(Re(eigen(Projection_matrix)$vectors[,1]))
right_eigenvector = right_eigenvector / sum(right_eigenvector)

classes = c('Seeds', 'Seedlings', 'Saplings', 'Poles', 'Mature trees')

predictions = 
  tibble(
    class = classes,
    population = right_eigenvector
  ) 

plot_stable_class_distribution =
  predictions |>
  ggplot(aes(class, population)) +
  geom_col(fill = c('salmon', "#DD5144")[1 + 1:5 %% 2]) +
  coord_flip() +
  theme(aspect.ratio = 1) +
  labs(
    y = 'Population (millions)',
    x = 'Age group'
  ) +
  ggtitle('Stable class distribution')


Reproductive value

As we discussed in lecture, ecologists / ecosystem managers / conservationists want to know not only the stable class distribution and growth rate of the population, but also how much each class contributes to the growth of every other class. In other words, we want to know the reproductive value of each class. For example, conservationists may use this information to determine which class is the most critical to protect for the health of the population.

Just as the long-term relative abundances of each class are proportional to the components of the leading right eigenvector of the projection matrix, the reproductive value of each class is proportional to the components of the leading left eigenvector.

Note: left and right eigenvectors:

  • A right eigenvector of a matrix \(\textbf{M}\) is by definition a vector \(\textbf{v}\) that satisfies the equation \(\textbf{M}\textbf{v} = \lambda \textbf{v}\), where \(\lambda\) is an eigenvalue of \(\textbf{M}\).

  • A left eigenvector of a matrix \(\textbf{M}\) is by definition a vector \(\textbf{u}\) that satisfies the equation \(\textbf{M}^\textrm{T} \textbf{u}= \lambda \textbf{u}\), where \(\textbf{M}^\textrm{T}\) is the transpose of \(\textbf{M}\).

Note: matrix transpose:

  • The transpose of a matrix \(\textbf{M}\) is the matrix \(\textbf{M}^\textrm{T}\) defined such that \(M^\textrm{T}_{ij} = M_{ji}\). In words, we flip the matrix so its rows become columns and vice-versa.

  • The eigenvalues of \(\textbf{M}^\textrm{T}\) are the same as the eigenvalues of \(\textbf{M}\).


The code chunk below computes the leading left eigenvector of the projection matrix above, then rescales it so its components add up to 1:

left_eigenvector = abs(Re(eigen(t(Projection_matrix))$vectors[, 1]))
left_eigenvector = left_eigenvector / sum(left_eigenvector)


We can now add this to our predictions tibble and write the plotting code.


Q10. Use the code below to plot the reproductive value of the different classes. Compare your results with the stable class distribution in Q9.

predictions %<>%
  mutate(reproductive_value = left_eigenvector)

plot_reproductive_value =
  predictions |>
  ggplot(aes(class, reproductive_value)) +
  geom_col(fill = c('deepskyblue1', "#4C8BF5")[1 + 1:5 %% 2]) +
  coord_flip() +
  theme(aspect.ratio = 1) +
  labs(
    y = 'Reproductive value',
    x = 'Age group'
  ) +
  ggtitle('Reproductive value')


Elasticity

Finally, another important bit of information is how the population would respond to a change in its vital rates. For example, a fisheries manager needs to know what percentage of the adult population they can harvest (which reduces adult survival) without jeopardizing the sustainability of the system. Or an evolutionary biologist may be interested in how a change in the reproductive schedule of a species (ie changing fecundity rates of the different classes) would affect its fitness (ie population growth rate).

Several statistics can been used to answer these questions. Here we will look at the elasticity of the system, a concept borrowed from economics which measures the percentage change in the population growth upon a percentage change in each of the vital rates.

The mathematical details need not concern us, but this paragraph will help understand the code below. The elasticity is formally defined as \(e_{ij} = \frac{\partial \log(\lambda)}{\partial \log(P_{ij})}\). In words, it is the rate of change of the logarithm of the leading eigenvalue upon a small increment in the logarithm of the entries of the projection matrix. The good news is that it can be directly obtained from quantities we have already looked at, namely the matrix entries \(P_{ij}\), the leading eigenvalue \(\lambda\), and the left and right leading eigenvectors \(l_i, \; r_i\), via the formula \(e_{ij} = \frac{1}{\lambda}P_{ij} \, l_i \, r_j\).

# formula for elasticity
elasticity = 
  outer(left_eigenvector, right_eigenvector) * 
  Projection_matrix / 
  leading_eigenvalue

# renormalize so elements add up to 1
elasticity = elasticity / sum(elasticity)
  
# assemble a tibble for plotting
elasticity_tibble = 
  tibble(
    Pij = as.numeric(Projection_matrix),
    elasticity = as.numeric(elasticity),
    row = factor(rep(classes, times = length(classes)), levels = rev(classes)),
    column = factor(rep(classes, each = length(classes)), levels = classes)
  )

# plotting code
plot_elasticity = 
  elasticity_tibble |>
  ggplot(aes(column, row, fill = elasticity)) +
  geom_raster() +
  scale_x_discrete(position = "top") +
  geom_text(
    aes(
      x = rep(1:5, each = 5), 
      y = rep(5:1, times = 5)
    ), 
    label = as.numeric(round(elasticity, 2)), 
    color = 'white',
    fontface = 'bold'
  ) +
  theme(
    aspect.ratio = 1,
    legend.position = 'none'
  ) +
  ggtitle('Elasticity')

plot_elasticity

The tiles above show how strongly the population growth rate responds to an increment in the corresponding entries of the projection matrix. We conclude that increasing the survival of mature trees has by far the strongest impact on population growth, followed by increasing survival of poles and saplings. By comparison, increasing fecundity of mature trees has minimal impact. This was not at all obvious from the projection matrix, especially since fecundity and survival have different units (seeds per individual and proportion of surviving individuals), and hence are not directly comparable.

This is precious information if we wish to protect American beeches. If we had to decide, for example, between spending our resources on

  • Option 1: Increasing mature tree survival by a certain percentage – e.g. by reducing the prevalence of beech bark disease via controls on the population of the beech scale insect, or

  • Option 2: Increasing mature tree fecundity by the same percentage – e.g. by protecting beechnuts from seed predators,

then our results above indicate that Option 1 is almost 30 times more effective at increasing population growth of this species.


Figure. Cryptococcus fagisuga (left), the beech scale insect, damages the bole of F. grandifolia by boring holes in its bark to feed on sap, leaving wounds that are later infected by the beech bark disease fungus Neonectria faginata (right), which can ultimately kill the tree. Pictured are insects covered in a white wool-like protective secretion. Photos by John A. McLaughlin.


Conversely, suppose we are evaluating an American beech stand for potential commercial exploitation. We can either harvest mature trees for timber, or seeds for the food industry (beechnuts are edible but too small to be commercially viable, but let’s go with the hypothetical here). Based on the very high survival of mature trees, we might assume that harvesting 5% of the mature trees every year would have minimal impact on the health of the population. Alternatively, given the large seed production per tree, we might assume that we can safely harvest 5% of the nuts every year.


Q11. What is the projected impact of each alternative plan on the tree population? Run the code below, which simulates the population of mature trees per acre under each scenario above for 100 years, and describe the outcome. How do your results mesh with the elasticity results discussed above?

# function for time-evolving the population given a projection matrix 
# and initial abundances
Project = function(N0, P, years){
  N = NULL
  for(year in years){
    N = cbind(N, matrix.power(P, year) %*% N0)
  }
  return(N)
}

# setting up the different management scenarios
P_harvest_seeds = P_harvest_adults = P_unmanaged = Projection_matrix

# 5% harvest
P_harvest_adults[5, 5] = Projection_matrix[5, 5] * .95
P_harvest_seeds[1, 5] = Projection_matrix[1, 5] * .95

# initial abundances per acre of each class
N0 = c(20000, 2000, 300, 100, 100)

# time-evolve the population under each scenario, collect results under one tibble
N =
  tibble(
    year = 0:100,
    `Unmanaged` = Project(N0 = N0, P = P_unmanaged, years = year)[5, ],
    `Harvest 5% seeds` = Project(N0 = N0, P = P_harvest_seeds, years = year)[5, ],
    `Harvest 5% adults` = Project(N0 = N0, P = P_harvest_adults, years = year)[5, ]
  )

# plotting code
plot_N =
  N |>
  pivot_longer(-year, names_to = 'treatment', values_to = 'adults') |>
  ggplot(aes(year, adults, group = treatment, color = treatment)) +
  geom_line(size = 1) +
  theme(
    aspect.ratio = 1
  ) +
  ylab('Mature trees per acre')

# call the plot
plot_N


Q12. Through trial and error, estimate what percentage of seeds can be safely harvested from this population while keeping it from eventual collapse. You can either use the plots or the leading eigenvalue of your projection matrix to arrive at your estimate.


Post-hurricane recovery of Vochysia ferruginea

Figure: The showy yellow flowers of Vochysia ferruginea are a striking feature of the rainforest landscape. Photo by Reinaldo Aguilar

In the aftermath of Hurricane Joan in southeastern Nicaragua in 1988, researchers were concerned about the long-term viability of the local population of Vochysia ferruginea, a tree native to lowland tropical rainforests of South America. Having monitored the population dynamics of the species for five years following the Category 4 hurricane, Boucher and Mallona (1997) constructed the following projection matrix:

Table: Stage-projection matrix for V. ferruginea, from Boucher and Mallona (1997) via Vandermeer and Goldberg (2013)
Stage Seedling Small sapling Large sapling Small adult Large adult
Seedling 0.209 0.000 0.000 35.600 70.100
Small sapling 0.010 0.653 0.020 0.000 0.000
Large sapling 0.000 0.170 0.407 0.000 0.000
Small adult 0.000 0.000 0.570 0.731 0.000
Large adult 0.000 0.000 0.000 0.266 0.997


Q13. What is the biological meaning of the non-zero value of the entry on Row 2 (Small sapling), Column 3 (Large sapling) of the projection matrix? What phenomena may be responsible for this?


To answer questions Q14-Q17, you will need to find the leading eigenvalue and the leading right and left eigenvectors of the projection matrix P below. You can do that using the following lines of code:

leading_eigenvalue = Re(eigen(P)$values[1])
right_eigenvector = abs(Re(eigen(P)$vectors[, 1]))
left_eigenvector = abs(Re(eigen(t(P))$vectors[, 1]))


Q14. Based on this projection matrix (code below), determine whether this population is projected to grow or decline.

classes = c('Seedling', 'Small sapling', 'Large sapling', 'Small adult', 'Large adult')

P = 
  matrix(
    c(
      .209,    0,    0, 35.6, 70.1,
       .01, .653,  .02,    0,    0,
         0,  .17, .407,    0,    0,
         0,    0,  .57, .731,    0,
         0,    0,    0, .266, .997
    ),
    byrow = TRUE,
    nrow = 5
  )


Q15. Between small adults and large adults, which class is projected to be more abundant in the long term?


Q16. Which class has the highest reproductive value?


Q17. Plot the elasticity matrix for this population. Identify which of the vital rates in the projection matrix would lead to the biggest change in population growth if modified. Hint: You will need to copy/paste and modify the code chunk given above that plots the elasticity matrix for American beech.



—

Appendix: Matrix algebra, Eigenvalues, and Eigenvectors

Class-structured population growth as a \(2^+\)-dimensional analogue of exponential growth

This appendix will help you make intuitive sense of vectors, matrices, eigenvalues, and eigenvectors in case you have little or no previous exposure to them.

Consider the simplest type of population growth model, exponential growth:

\[ N_{t+1} = \lambda N_t \]

The per capita growth rate \(\lambda\) is a composite of the birth rate \(b\) (also called fertility or fecundity) and the death rate \(d\) (aka mortality). Specifically, we can break down the expression above into its component processes:

\[ \begin{eqnarray} N_{t+1} &=& N_t + b N_t - d N_t \\ &=& (1-d)N_t + b N_t \\ &=& pN_t + bN_t \end{eqnarray} \]

The first line says that the population next year (or any other time unit) is equal to its current size this year plus the new births minus the deaths, and the second line gathers the first and last term, and says next year’s population will be those who survived from this year plus the new arrivals. In the last line we simply define \(p = 1-d\) as the survival rate of the population from one year to the next. In the end, we have a very simple relationship between the net growth rate of the population \(\lambda\), which is a population-level quantity, and the individual-level rates of survival \(p\) and reproduction \(b\), \(\lambda = p + b\).

This week, we introduced one layer of complexity, namely class structure in the population. We abandon the simplistic construct of one-size-fits-all birth and death rates and recognize that

  • Each class \(k\) has its own fecundity \(b_k\) and its own mortality \(d_k\) (or equivalently its own survival),

  • Births can only increase the population of the first class (since every individual begins life in the first class),

  • Individuals may transition between classes over time.

Consider a population containing two classes, say juveniles \(J\) and adults \(A\). Assuming all birth, death, and transition rates depend only on class and not on the environment or the current size of each class (which is an unrealistic assumption in many circumstances, as we will see later in the course), our model now takes the form

\[ J_{t+1} = p_J J_t + b_J A_t \\ A_{t+1} = a_J J_t + p_A A_t \tag{1} \]

In words: the number of juveniles next year will be the number of juveniles who survived this year plus the number of new births (which is proportional to the number of egg-laying or birth-giving adults this year), while the number of adults next year is the number of adults who survived this year plus the number of juveniles who transitioned into adulthood (which is proportional to the number of juveniles this year).

This may look a lot more complicated than \(N_{t+1} = (p+b) N_t\), but really it is saying the same thing: we have a survival term and a new arrivals term, except that now rather than keep track of a single population, we are now keeping track of two connected populations.

In particular, both equations above have a term proportional to \(J\) and a term proportional to \(A\). This is very fortunate, because there are handy mathematical shortcuts we can take with systems of equations of this form. (It may not seem like it, but being able to borrow from well-established fields of mathematics simplifies things enormously! This is the reason why we talk about matrices and eigenvalues and eigenvectors in a biology class.)

We can write the equations above in matrix format:

\[ \pmatrix{J_{t+1} \\ A_{t+1}} = \pmatrix{p_J & b_A \\ a_J & p_A} \pmatrix{J_{t} \\ A_{t}} \tag{2} \]

where we grouped the populations of juveniles and adults into a two-component vector \(\pmatrix{J_{t} \\ A_{t}}\), which you can think of as the two-dimensional extension of the one-dimensional population size \(N\). If we introduce the vector notation \(\bf{N}_t \equiv \pmatrix{J_{t} \\ A_{t}}\) (notice the boldface \(\bf{N}\) to distinguish it from a regular one-dimensional number \(N\)), and represent the projection matrix \(\pmatrix{p_J & b_A \\ a_J & p_A}\) by the symbol \(\bf{P}\), then the analogy with the unstructured model of population growth becomes clear:

\[ \bf{N}_{t+1} = P \bf{N}_t \]

This expression looks very much like (one-dimensional) exponential growth \(N_{t+1} = \lambda N_t\): the (two-dimensional) population next year is equal to its current size multiplied by a (group of) constant(s), which must work as some kind of growth rate. So just as \(\bf{N}_t\) is the 2-dimensional extension of the population, the projection matrix \(\bf{P}\) is the 2-dimensional extension of the growth rate \(\lambda\).

To make precise this intuitive analogy between two-dimensional class-structured population growth and one-dimensional exponential growth, we must learn a bit of matrix algebra:

A casual primer on matrix algebra

When you multiply a two-dimensional number—ie a vector— \(\bf{N} = (\it{J}, \it{A})\), by a constant \(\lambda\), you get another vector, \(\lambda\bf{N} = (\lambda \it{J}, \lambda \it{A})\). Notice that each component of the new vector is proportional to the original component via the same proportionality constant, \(\lambda\). For those familiar with the arrow representation of vectors, this means that the new vector points in the same direction as the old one.

When you multiply a vector by a matrix, you also get another vector, just like when you multiply the vector by a constant. However, the difference is that typically the new vector does not point in the same direction as the old one; in other words, the new components are not simply proportional to the old ones. Instead, you get a mix of the old components multiplied by different constants.

The formal rules for multiplying a vector by a matrix are

\[ \pmatrix{a_{11} & a_{12} \\ a_{21} & a_{22}} \pmatrix{x \\ y} = \pmatrix{a_{11} x + a_{12} y \\ a_{21} x + a_{22} y} \]

Notice that our two-dimensional (ie two-class) structured population growth model follows this rule, as you can verify by examining equations (1) and (2) above. Biologically, our observation that the new vector is not proportional to the old one corresponds to observing that the proportion of juveniles to adults changes from one year to the next.

Why do we care about this? Let’s assume for a second that it was indeed the case that the new populations of juveniles and adults were simply equal to the old ones multiplied by a constant, say, \(\lambda\). Mathematically, we would write \(\bf{N}_{t+1} = P \bf{N}_t = \lambda N_t\), which in components notation would be written \(J_{t+1} = \lambda J_t\), \(A_{t+1} = \lambda A_t\). That would mean that the ratio between the number of juveniles and adults remains unchanged from one year to the next, \(J_{t+1} / A_{t+1} = J_t / A_t\). So if in Year \(t\) we had a 2-to-1 juveniles-to-adults ratio, in Year \(t+1\) we would still have a 2-to-1 juveniles-to-adults ratio. The only difference is that the population as a whole will have grown by a factor \(\lambda\) which is of course the growth rate of the population!

What would happen in Year \(t+2\)? We would have \(\bf{N}_{t+2} = P \bf{N}_{t+1} = P(\lambda N_t) = \lambda (P N_t) = \lambda\lambda N_t = \lambda^2 N_t\). In other words, the juveniles-to-adults ratio remains constant in Year \(t+2\) as well, except by now the population has grown by a factor \(\lambda^2\). Extending this logic forward, we see that population structure (ratio of juveniles to adults) would be permanent, and the population as a whole will be growing exponentially with annual rate \(\lambda\) – just like in the much simpler one-dimensional model! This permanent population structure is called the stable age/stage/class distribution.

Interestingly, it turns out that, given a projection matrix \(\bf{P}\), we can find a special vector \(\bf{v}\) whose components are in such proportion that the vector \(\bf{Pv}\) is proportional to \(\bf{v}\). Let’s look at an example. Consider the projection matrix

\[ \bf{P} = \pmatrix{0.1 & 10 \\ 0.84 & 0.9} \]

In a population with this projection matrix, 10% of juveniles survive till next year, while 84% transition to adulthood, and meanwhile each adult produces on average 10 offspring and have a 90% of surviving till next year.

Suppose the population currently has 300 juveniles and 100 adults. What will this population look like next year? We find that out by multiplying the vector \(\pmatrix{300 \\ 100}\) by the projection matrix:

\[\pmatrix{0.1 & 10 \\ 0.84 & 0.9}\pmatrix{300 \\ 100} = \pmatrix{1030 \\ 343} = \mathrm{3.43}\pmatrix{300 \\ 100}\]

This is certainly not true of every initial state. For example,

\[\pmatrix{0.1 & 10 \\ 0.84 & 0.9}\pmatrix{200 \\ 100} = \pmatrix{1020 \\ 258} = 2.58\pmatrix{395 \\ 100}\]

In concrete terms, this means that if a population with this projection matrix has a 3-to-1 juveniles to adults ratio this year, then it will also have a 3-to-1 ratio next year, and every year after. But if the current ratio is not 3, then it will not remain the same next year, so the class structure of the population will be changing from one year to the next.

Because the vector \(\pmatrix{300 \\ 100}\) has this special relationship with the matrix , we call it an eigenvector of that matrix. And because multiplying the matrix by that vector is equivalent to simply multiplying the vector by \(3.43\), we say that \(3.43\) is the eigenvalue corresponding to that eigenvector (from the German eigen, meaning “own”). Notice that the eigenvalue ends up being the growth rate of the population.

In general, a two-dimensional projection matrix like the one above will have not one but two eigenvectors, each with its own eigenvalue. In the case of our example matrix, the other eigenvector is \(\pmatrix{-396 \\ 100}\) and the corresponding eigenvalue is \(-2.42\) (don’t worry about the negative signs). More generally, an \(n\)-dimensional matrix has \(n\) eigenvectors and \(n\) eigenvalues.

By now you can probably see why (one of the) eigenvector(s) of the projection matrix predicts the stable class distribution (in our example, the long-term juveniles-to-adults ratio). But assuming the population is not currently in its stable class distribution, what reason do we have to believe that it will ever reach it?

The reason is another fortunate mathematical fact: every two-dimensional vector is a mixture of the two eigenvectors of the projection matrix. In matrix notation, we can write

\[ \pmatrix{J \\ A} = c_1 \pmatrix{300 \\ 100} + c_2 \pmatrix{-396 \\ 100} \] where the coefficients \(c_1\) and \(c_2\) make the equality exact. For example, for the vector \((J=200,A=100)\), which as we have seen is not an eigenvector of \(\bf{P}\), we have \(c_1 = 0.856\) and \(c_2 = 0.144\). The coefficients can be interpreted as weights: the vector \((J=200, A=100)\) is 85.6% proportional to the eigenvector \((J=300, A=100)\) and 14.4% proportional to the eigenvector \((J=-396,A=100)\), or approximately 6 times more similar to the first eigenvector than the second. (Notice here I am using a more inline-friendly vector notation, which in the next paragraph I’ll streamline further by dropping the letters \(J\) and \(A\).)

What happens if we multiply \((200, \, 100)\) by our projection matrix (ie if we advance a starting population with 2 juveniles and 1 adult to the next year)? As seen above, we get \((1020, \, 258)\). How does this new vector relate to the eigenvectors? Some algebra reveals that the coefficients here are \(c_1 = 2.93\) and \(c_2 = -0.35\). Ignoring the negative sign, this means that now the vector is 89.2% proportional to \((300, \, 100)\) and only 10.8% proportional to \((-396, \, 100)\), or approximately 8.3 times more similar to the first eigenvector than the second. Another way to put this is that after one year the juveniles-to-adults ratio is now more similar to the 3-to-1 ratio of the first eigenvector than it was before.

If you advance the population one more year, ie multiply this resulting vector by the projection matrix again, the population becomes almost 12 times more similar to the first eigenvector than the second. As the years go by, that similarity ratio hikes up to 16, 24, 32, and so on forever upwards. Clearly, the population comes to look more and more like the first eigenvector. And that is not a coincidence! The reason the first eigenvalue eventually dominates is because its eigenvalue has a larger magnitude than the eigenvalue of the other eigenvector. To see this, let’s represent the first and second eigenvectors by the symbols \(\bf{v}_1\) and \(\bf{v}_2\), and their respective eigenvalues as \(\lambda_1\) and \(\lambda_2\), noting that \(|\lambda_1| > |\lambda_2|\). In Year 0, our population is

\[ \bf{N}_0 = \mathrm{c}_1 \bf{v}_1 + \mathrm{c}_2 \bf{v}_2 \] That is, in Year 0 the similarity of the population structure to \(\bf{v}_1\) relative to its similarity to \(\bf{v}_2\) is \(\frac{|c_1|}{|c_2|}\).

In Year 1, the population will be

\[ \begin{eqnarray} \bf{N}_1 &=& \bf{PN}_0 \\ &=& \bf{P}(\mathrm{c}_1 \bf{v}_1 + \mathrm{c}_2 \bf{v}_2) \\ &=& \mathrm{c}_1 \, \bf{Pv}_1 + \mathrm{c}_2 \, Pv_2 \\ &=& \mathrm{c}_1 \lambda_1 \bf{v}_1 + \mathrm{c}_2 \lambda_2 \bf{v}_2 \end{eqnarray} \]

which means that in Year 1 that similarity proportion will be \(\frac{|\lambda_1|}{|\lambda_2|}\frac{|c_1|}{|c_2|}\). Since \(|\lambda_1| > |\lambda_2|\), we can see that the population becomes more similar to \(\bf{v}_1\), regardless of the values of \(c_1\) and \(c_2\).

In Year 2, the population will be

\[ \begin{eqnarray} \bf{N}_2 &=& \bf{PN}_1 \\ &=& \bf{P}(\mathrm{c}_1 \lambda_1 \bf{v}_1 + \mathrm{c}_2 \lambda_2 \bf{v}_2) \\ &=& \mathrm{c}_1 \lambda_1 \, \bf{Pv}_1 + \mathrm{c}_2 \lambda_2 \, Pv_2 \\ &=& \mathrm{c}_1 \lambda_1\lambda_1 \bf{v}_1 + \mathrm{c}_2 \lambda_2\lambda_2 \bf{v}_2 \\ &=& \mathrm{c}_1 \lambda_1^2 \bf{v}_1 + \mathrm{c}_2 \lambda_2^2 \bf{v}_2 \end{eqnarray} \]

In other words, in Year 2 the population will be \(\frac{|\lambda_1|^2}{|\lambda_2|^2}\frac{|c_1|}{|c_2|}\) times more similar to \(\bf{v}_1\) than \(\bf{v}_2\).

By now you can probably see the pattern. In Year \(t\), the population will be \(\bf{N}_t = \mathrm{c}_1 \lambda_1^t \bf{v}_1 + \mathrm{c}_2 \lambda_2^t \bf{v}_2\), and the similarity proportion will be \(\frac{|\lambda_1|^t}{|\lambda_2|^t}\frac{|c_1|}{|c_2|} = \left(\frac{|\lambda_1|}{|\lambda_2|}\right)^t\frac{|c_1|}{|c_2|}\). And since the ratio \(\frac{|\lambda_1|}{|\lambda_2|}\) is greater than 1, it only increases when raised to a power. The higher the power, the closer the similarity to the first eigenvector.

We conclude that no matter what the initial stage of the population, over time it will look more and more like the eigenvector whose eigenvalue has the higher magnitude. This larger eigenvalue is called the leading eigenvalue, and the corresponding eigenvector is the leading eigenvector.

Now you know why the leading eigenvector predicts the stable class distribution of the structured population.

NOTE: Everything we discussed above for a two-class population and its corresponding two-dimensional projection matrix, the two different eigenvectors, etc, extends to higher number of classes in the expected way.


We close this discussion by pointing out another special vector associated with a projection matrix and its biological meaning. Note that the entry on row \(i\) and column \(j\) of the projection matrix tells us what contribution the next cohort of Class \(i\) receives from the current cohort of Class \(j\). It is also convenient to introduce the transpose of the projection matrix, symbolized \(\bf{P}^T\), where we flip the rows and columns of \(\bf{P}\): \(p^T_{ij}=p_{ji}\). For our example above, we have \(\bf{P}^T=\pmatrix{0.1 & 0.84 \\ 10 & 0.9}\). Notice that the element on row \(i\) and column \(j\) of the transpose projection matrix tells us how Class \(i\) contributes to Class \(j\). The transpose matrix \(\bf{P}^T\) will have its own eigenvectors. The following are some facts about those eigenvectors, which I will state without proof (for a less cursory discussion of this, see Vandermeer and Goldberg 2013, and for an in-depth treatment see Caswell 2001).

  • A matrix has right eigenvectors and left eigenvectors. The eigenvectors we have been discussing in this Appendix are right eigenvectors. When the word eigenvector is used without qualification, it is implied that one is speaking of right eigenvectors.

  • If \(\bf{v}\) is a right eigenvector of a matrix \(\bf{P}\), then \(\bf{Pv} = \lambda\bf{v}\).

  • If \(\bf{u}\) is a left eigenvector of a matrix \(\bf{P}\), then \(\bf{P^Tu} = \lambda\bf{u}\).

  • In other words, the right eigevenctors of \(\bf{P}^T\) are the left eigenvectors of \(\bf{P}\).

The reason we are bothering with this other type of eigenvector is that the leading left eigenvector of \(\bf{P}\) describes the reproductive value of each class. This concept is a holistic summary of the fertility of the class plus its expected fertility later when the individual may have transitioned to other classes. In short, it measures how much the class is expected to contribute to the growth of the population.